program main

    use gauss
    use exp, only: new_solve => solve, solve_f => solve1
    real(kind=8), allocatable :: A(:, :), b(:), c_(:), c1(:), b_(:,:)
    real :: tic, toc
    
#ifdef HUGE
    integer, parameter :: N = 2500
#else
    integer, parameter :: N = 500
#endif
    integer :: ipiv(N)

    allocate(A(N, N), b(N), c_(N), c1(N), b_(N,1))

    call random_number(A)
    call random_number(b)

    ! call cpu_time(tic)
    ! c_ = linear_solve(A, b)
    ! call cpu_time(toc)
    ! print "(A, ES12.4)", "solution1(UserDefined1), toc - tic(seconds): ", toc - tic

    ! call cpu_time(tic)
    ! c1 = linear_solve1(A, b)
    ! call cpu_time(toc)
    ! print "(A, ES12.4)", "solution2(UserDefined2), toc - tic(seconds): ", toc - tic

    ! call cpu_time(tic)
    ! call solve(A, b, N, 1)
    ! call cpu_time(toc)
    ! print "(A, ES12.4)", "solution3(Euler)       , toc - tic(seconds): ", toc - tic
    ! ! print "(A, ES12.4)", "L_2 norm: ", norm2(c1 - c_)

    b_(:, 1) = b
    call cpu_time(tic)
    call DGESV( N, 1, A, N, IPIV, b_, N, INFO )
    call cpu_time(toc)
    print "(A, ES12.4)", "solution4(LAPACK)      , toc - tic(seconds): ", toc - tic
    ! print "(A, ES12.4)", "L_2 norm: ", norm2(c1 - c_)

    ! call cpu_time(tic)
    ! call new_solve(A, b_)
    ! call cpu_time(toc)
    ! print "(A, ES12.4)", "solution5(Euler1)      , toc - tic(seconds): ", toc - tic
    
    call cpu_time(tic)
    b_ = solve_f(A, b_)
    call cpu_time(toc)
    print "(A, ES12.4)", "solution6(Euler2)      , toc - tic(seconds): ", toc - tic

end program main
